function [ L, D, U, P ] = ldu_decomp (A)
[m, n] = size(A);
p = 1:m;
for i=1:n
    A(i:m, i)
    [mv, mi] = max(abs(A(i:m, i)));
%     mi = mi + i - 1;
%     if mi ~= i
%         display('pw')
%         %Need to make a permutation
%         %Swap the row
%         rowtoswap = A(i , :);
%         A(i , :) = A(mi , :);
%         A(mi, :) = rowtoswap;
%         %Swap the index
%         previouspivot = p(i);
%         p(i) = p(mi);
%         p(mi) = previouspivot;
%     end
end
L = eye(m,m);
D = eye(m,m);
U = A;
for i=1:n
    for j=(i+1):n
        if U(i, i) ~=0
            L(j, i) = U(j, i) ./ U(i, i);
        else
            L(j, i) = 0;
        end
        U(j, i:n) = U(j, i:n) - L(j, i) * U(i ,i:n);
    end
    D(i, i) = U(i, i);
    if D(i , i) ~=0
        U(i, i:n) = U(i, i:n) ./ D(i, i);
    else
        U(i, i:n) = 0;
    end
end
% Create P from indices
P = zeros(m, n);
for i=1:m
    P(i, p(i)) = 1;
end